capture log close
log using ncreadexample.log, replace

* Example1 :Setting up the Java dependence
netcdf_init, download plus(jar)

* Example 2: Display the Metadata with ncread
////The developed commands can directly read nc files on the network. However, due to reasons such as network SSL authentication, the reading may fail. If this happens, you can copy the nc file to the local device and then perform the following corresponding operations.
local url = "https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/" + ///
"NEX-GDDP-CMIP6/BCC-CSM2-MR/ssp245/r1i1p1f1/tas/" + ///
"tas_day_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_2050.nc"
ncread using `"`url'"'

* Example 3: Display variable metadata with ncdisp
///tas variable
local url = "https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/" + ///
"NEX-GDDP-CMIP6/BCC-CSM2-MR/ssp245/r1i1p1f1/tas/" + ///
"tas_day_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_2050.nc"
ncdisp tas using `"`url'"'

///time variable
local url = "https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/" + ///
"NEX-GDDP-CMIP6/BCC-CSM2-MR/ssp245/r1i1p1f1/tas/" + ///
"tas_day_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_2050.nc"
ncdisp time using `"`url'"'

* Example 4: import grid data into Stata
///load grid data on 1st Jan 2050
local url = "https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/" + ///
"NEX-GDDP-CMIP6/BCC-CSM2-MR/ssp245/r1i1p1f1/tas/" + ///
"tas_day_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_2050.nc"

ncread tas using `"`url'"' , origin(1 1 1) size(1 -1 -1) clear

sum

list in 1/10

///load grid data that encompasses Hunan Provinc
local url = "https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/" + ///
"NEX-GDDP-CMIP6/BCC-CSM2-MR/ssp245/r1i1p1f1/tas/" + ///
"tas_day_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_2050.nc"

ncread lon using `"`url'"', clear 
gen n=_n 
qui sum n if lon>=108 & lon<=115
local lon_start = r(min)
local lon_count = r(N)

ncread lat using `"`url'"', clear 
gen n=_n 
qui sum n if lat>=24 & lat<=31
local lat_start = r(min)
local lat_count = r(N)

ncread tas using `"`url'"', clear origin(1 `lat_start' `lon_start') ///
 size(-1 `lat_count' `lon_count')
 
gen date = time - 3650.5  + date("2050-01-01", "YMD")
format date %td

list in 1/10

save "grid_all.dta", replace

* Example 5:Match cities to nearest four grid cells using matchgeop
use "grid_all.dta", clear

rename lon ulon
rename lat ulat
save "grid_all_1.dta", replace

keep if time==3650.5
gen n=_n 
save "hunan_grid.dta", replace

use "hunan_city.dta", clear
matchgeop ORIG_FID lat lon using hunan_grid.dta, neighbors(n ulat ulon) nearcount(4) gen(distance) bearing(angle)

merge m:1 n using hunan_grid.dta, keep(3)
drop _merge
save "hunan_origin.dta", replace

drop time date
joinby ulat ulon using grid_all_1.dta

sort ORIG_FID date 
list city distance angle date tas in 1/10

///Diagram of azimuthal angle
use "hunan_origin.dta", clear
keep if city == "Changsha"

sort angle
gen id=_n

local R = 6371 
gen lat_rad = lat * (_pi/180)  

gen delta_lat = (distance / `R') * (180/_pi)  
gen delta_lon = (distance / (`R' * cos(lat_rad))) * (180/_pi)  

expand 90
bysort n: gen t = _n - 1
gen theta = (angle * t/89) * (_pi/180)  

gen arc_lat = lat + delta_lat * cos(theta)
gen arc_lon = lon + delta_lon * sin(theta)

bysort n: gen label_theta = (angle/2) * (_pi/180)
gen label_lat = lat + delta_lat * cos(label_theta)
gen label_lon = lon + delta_lon * sin(label_theta)
replace label_lat = lat + delta_lat * 0.7 * cos(label_theta) if id == 3
replace label_lon = lon + delta_lon * 0.7 * sin(label_theta) if id == 3

gen latlon_label = "(" + string(lat, "%8.2f") + "°N, " + string(lon, "%8.2f") + "°E)"
gen ulatlon_label = "(" + string(ulat, "%8.2f") + "°N, " + string(ulon, "%8.2f") + "°E)"
gen angle_label = string(angle, "%8.2f") + "{&degree}"

twoway pcarrowi 28.15 113.15307 28.52 113.15307, color(black) ///
    || pcarrow lat lon ulat ulon, color(blue) ///
    || scatter lat lon if t == 1, mcolor(red) mlabel(latlon_label) mlabcolor(black) mlabpos(9) mlabgap(0.8) mlabsize(medium) ///
    || scatter ulat ulon if t == 1 & (ulon == 113.125), mcolor(red) mlabel(ulatlon_label) mlabcolor(black) mlabpos(9) mlabgap(0.8) mlabsize(medium) ///
    || scatter ulat ulon if t == 1 & (ulon == 113.375), mcolor(red) mlabel(ulatlon_label) mlabcolor(black) mlabpos(3) mlabgap(0.8) mlabsize(medium) ///
    || line arc_lat arc_lon if id == 1 , lcolor(green) ///
    || line arc_lat arc_lon if id == 2, lcolor(green) ///
    || line arc_lat arc_lon if id == 3, lcolor(green) ///
    || line arc_lat arc_lon if id == 4, lcolor(green) ///
    || scatter label_lat label_lon, mlabel(angle_label) msymbol(i) mlabcolor(black) mlabsize(medium) ///
    xscale(off noline) yscale(off noline) xlabel(, nogrid noticks) ylabel(, nogrid noticks) ///
    aspect(1) legend(off)

graph save gragh1, replace
	
* Example 6 : Calculate 80km-radius IDW temperatures for cities
use "grid_all.dta", clear
rename lon ulon
rename lat ulat
gen n=_n
save "grid_all_2.dta", replace

use "hunan_city.dta", clear
matchgeop ORIG_FID lat lon using grid_all_2.dta, neighbors(n ulat ulon) within(80) gen(distance)

merge m:1 n using grid_all_2.dta, keep(3)
drop _merge

list city ulat ulon distance date tas in 1/10

save "hunan_80km.dta", replace

drop if tas==.
gen weight  = 1/distance
gen weighted_tas = tas * weight
bysort city date : egen sum_weighted_tas = total(weighted_tas)
bysort city date : egen total_weight = total(weight)
gen idw_tas = sum_weighted_tas / total_weight
gen temp_c = idw_tas - 273.15

duplicates drop city date , force

list city  distance date  temp_c in 1/10

save "hunan_IDW.dta", replace

///Distribution of raster points within 80 kilometers of Changsha City
use "hunan_80km.dta" ,clear
keep if time == 3650.5
keep if city == "Changsha"

summarize lon 
local clon = r(mean)
summarize lat 
local clat = r(mean)

preserve
clear
set obs 720
gen theta = (_n-1) * 2 * _pi/720 

local R = 6371  
local d = 80    
local delta = `d'/`R'  

gen lat_c = asin(sin(`clat'*_pi/180)*cos(`delta') + ///
             cos(`clat'*_pi/180)*sin(`delta')*cos(theta)) * (180/_pi)
             
gen lon_c = `clon' + ///
            atan2(sin(theta)*sin(`delta')*cos(`clat'*_pi/180), ///
                  cos(`delta') - sin(`clat'*_pi/180)*sin(lat_c*_pi/180)) * (180/_pi)
				  
expand 2 if _n == _N
replace theta = 0 if _n == _N

save circle.dta ,replace
restore

merge 1:1 _n using circle.dta, nogen

twoway (line lat_c lon_c, sort lcolor(dimgray) lwidth(0.6) lpattern(solid)) ///
    (pcarrowi  28.22691 113.1531 28.22691 113.9696, color(black)) ///
	(scatter lat lon, mcolor(red) msymbol(D) msize(small)) ///
    (scatter ulat ulon, mcolor(blue) msymbol(O) msize(small)), ///
    xlabel(minmax) ylabel(minmax) ///
	aspect(1)  legend(off) ///
    xscale(off noline) yscale(off noline) xlabel(, nogrid noticks) ylabel(, nogrid noticks) ///
	text(28.26 113.71 "80km", size(medsmall) justification(center))

graph save gragh2, replace

///IDW interpolated temperature distributions in Hunan
shp2dta using "hunan.shp", database(hunan_db) coordinates(hunan_coord) genid(id)
//January 1

use "hunan_IDW.dta" ,clear
keep if date == date("01jan2050", "DMY")

spmap temp_c using "hunan_coord.dta", id(OBJECTID) clmethod(q) cln(6) fcolor(Heat) legtitle("Temperature (°C)") title("Temperature(20500101)") subtitle("Within 80km Radius")  
graph save graph1, replace

//July 1
use "hunan_IDW.dta" ,clear
keep if date == date("01jul2050", "DMY")

spmap temp_c using "hunan_coord.dta", id(OBJECTID) clmethod(q) cln(6) fcolor(Heat) legtitle("Temperature (°C)") title("Temperature(20500701)") subtitle("Within 80km Radius")  
graph save graph2, replace

graph combine graph1.gph graph2.gph

graph save gragh3, replace


log close